On hydrodynamic flow of expanding Bose-Einstein condensates 
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We study expansion of quasi-one-dimensional Bose-Einstein condensate (BEC) after switching off 
the confining harmonic potential. Exact solution of dynamical equations is obtained in framework 
of the hydrodynamic approximation and it is compared with the direct numerical simulation of the 
full problem showing excellent agreement at realistic values of physical parameters. We analyze 
the maximum of the current density and estimate the velocity of expansion. The results of the ID 
analysis provides also qualitative understanding of some properties of BEC expansion observed in 
experiments. 
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PACS numbers: 

Many current experiments with Bose-Einstein conden- 
sates (BEC) include regime of a free expansion of the gas 
which initially was (magnetically or optically) trapped 
by a confining potential 0, 0] ■ The respective phase of 
the BEC evolution is of practical interest since much in- 
formation about coherent matter waves is experimentally 
obtained at this stage (say, by absorption imaging), and 
also because of importance of a free BEC flow in a con- 
text of recently proposed new devices like, for example, 
an atomic laser [3j. At ultra-low temperatures, a trapped 
BEC is well described by the three-dimensional (3D) 
Gross-Pitaevskii (CP) equation which in the absence of a 
trap potential is known as a nonlinear Schrodinger (NLS) 
equation. In a number of specific cases it can be reduced 
to a quasi-one-dimensional one. More specifically, this 
happens when one can neglect (due to some reasons) in- 
teraction between transverse modes. In the absence of 
a trap the problem thus reduces to an integrable model, 
ID NLS equation, which is very well studied fundamental 
model of the nonlinear physics. 

The problem of theoretical description of evolution of 
BEC confined by harmonic potential with varying pa- 
rameters has been addressed in several papers |^|. In 
the experiment |]| the realization of BEC expansion in 
quasi- ID waveguide has been reported and again excel- 
lent agreement has been found with theoretical predic- 
tions at long enough values of time of evolution. These re- 
sults have been confirmed by numerical solution of corre- 
sponding dynamical equations for quasi- ID case 
the present communication we give complete analytical 
treatment of this problem. We consider only BECs with 
positive scattering lengths. In terms of the NLS equation, 
the situation we are dealing with corresponds to the evo- 
lution of an initially localized pulse with the initial profile 
corresponding to the ground state of a BEC confined in 
a parabolic potential, i.e. to the defocusing NLS equa- 
tion with zero boundary conditions at infinity. As it is 
known, such problem does not have soliton solutions and 
independently of number of particles the condensate will 
spread out and its density will tend to zero with time. 
The formal analytical description of the respective solu- 
tion at t — ► oo was obtained rather long ago 0|. However, 



for practical purposes of analysis of experimental data it 
is desirable to have a description of the condensate evo- 
lution during the initial stages of time, too. It turns out 
that such evolution also admits rather complete analyt- 
ical description in a number of cases, and, in particular, 
when the initial density distribution is smooth enough, 
or in terms of BEC in the limit of large number of atoms 
which leads to the so-called Thomas-Fermi (TF) approx- 
imation. Then the hydrodynamic approach allows one to 
describe analytically evolution of the density and velocity 
fields of the condensate. The evolution of BEC at initial 
stages appears to be quite rich. For example, as it has 
been recently shown in Q , it may display wave breaking 
for some specific initial distributions of the BEC density. 

In the present report we describe the free expansion 
of a BEC initially confined in a harmonic potential. On 
the base of the hydrodynamic approach we find the time 
dependence of the density and velocity distributions and 
calculate such characteristics of the gas flow as the value 
and space and time coordinates of the maximum of cur- 
rent, velocity of the "edge" points of the condensate, 
as well as asymptotic density and velocity distributions. 
Also we show that simple analytical estimates of the ID 
problem allow one to give qualitative understanding of 
the phenomena observed in experiments with effectively 
2D BEC. 

Let us start with the 3D GP equation for the order 
parameter ip = VK r J^)i 

dib h 2 

ih-^ = -—k-iP + V trap {v)-il) + g Q \i>\ 2 ^, (1) 

where we use the standard notation: go = Airh 2 a s /m, a s 
being the s-wave scattering length, which is considered 
positive, m being the atomic mass; V trap (r) is a trap 
potential. Considering the case of a pancake BEC 
we take Vt ra p = ^w 2 x 2 where ui is the harmonic oscil- 
lator frequency and in the transverse direction (i.e. in 
the direction orthogonal to the x-axis) the size of the 
condensate is supposed large enough to be considered 
infinite in the first approximation. In longitudinal direc- 
tion the size of the condensate is of order of magnitude 
a = (h/muj) 1 / 2 . It is convenient to introduce some typ- 
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ical reference frequency loq so that the trap frequency is 
measured in units of ujq, uj — vlo q , v being the dimcn- 
sionless trap frequency. To make the dynamical equation 
dimensionless, we introduce 



V>(r,t) 
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exp ( ikj^r^ — i — —t] ^(x, t), 
2m J 

(2) 



where rj_ = (y, z) and a 2 , = h/moj^ = a 2 ^, and make 



a change of independent variables x — 2~ 1 / 4 clqx' , t = 
2 1 / 2 t' /luo. This results in the canonical form of the NLS 
equation with a parabolic potential 



a* a 2 * 



dt dx 2 



2|ilr| 2 «f = - V 2 X 2 ^, 



(3) 



where the primes were suppressed. The dimensionless 
BEC wave function ^(x,t) is normalized according to 



f 4- 
/ \y( x ,t)\ 2 dx = - 



Air Ma s aQ 



/4 5 



(4) 



where A/" is the total number of particles and S is an 
effective area of the transverse cross section of the con- 
densate. 

A stationary solution ty(x,t) — *&(x) exp(— ip£), corre- 
sponding to the ground state of BEC, is given by ^{x) 
satisfying the equation 



d 2 ^ 

dx 2 



- 2|*| 2 * = -v 2 x 2 ^, 



(5) 



subject to the zero boundary condition at \x\ — > co and 
having no other zeros. The eigenvalue p, (chemical po- 
tential) is determined by the normalization J3J. In di- 
mensionless units the longitudinal size of the condensate 
is of order of magnitude p 1 / 2 jv and if p v, the con- 
siderable part of the condensate can be described by TF 
approximation in which the term with the second space 
derivative in J5J can be neglected almost everywhere, so 
that 



Po(x) = \Vtf(x)\ 2 = \{p-\v 2 x 2 ) 



with normalization 



Zfi/v (2u) 3 / 2 
\^ TF {x)\ 2 dx- {m 

3/Z/v 



3v 



(6) 
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TF approximation © fails at the tails of the density dis- 
tribution where the density decays exponentially as \x\ — > 
oo instead of vanishing at finite distance \x\ — ^J2p/v 
(called TF radius) according to Eq. ©. Since only a 
small part of the condensate's mass is concentrated in 
these tails (for example, for p — 2 and v = 0.5, it is less 
then 1.5%). Thus in the limit p ^S> v the distribution © 
is assumed to be a good approximation of the initial den- 
sity distribution of BEC before switching off the external 
potential. 



By equating Eqs. (@J and Q we determine the value 
of the dimensionless chemical potential p in terms of ex- 
perimentally measurable parameters 



= 2 1 / 4 3t 
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Then the condition yields the following criterion 

Afa s a/S 3> 1 of applicability of the TF approximation in 
physical units. In what follows we deal mainly with two 
parameters v and p, which completely define the initial 
distribution of the condensate density. 

After switching off the potential, the condensate 
evolves according to Eq. with v — 0. If evolution 
does not lead to wave breaking of the pulse, then we 
still can neglect the dispersion effects and use the hydro- 
dynamic approximation for description of this evolution. 
As we shall see, this is the case of initial distribution JBJ, 
so that the hydrodynamic approximation is valid even 
at asymptotically large values of time considered in 0. 
Therefore we pass from the NLS equation J3J) to its hy- 
drodynamic representation. We represent ^>(x,t) in the 
form 



$(x, t) — y p(x, t) exp I i I v(x',t)dx' 



(9) 

so that substitution of © into © with v = yields 

\p t + (pv) x = 0, \v t + vv x + p x = 0, (10) 
which is subject to the the initial conditions 



p(x, 0) = p (x), v(x, 0) = 0, 



(11) 



where po(x) is taken to be the initial distribution JfjJ). In 
ffTTTfl we have neglected the higher space derivatives of the 
density p(x, t) that correspond to the quantum-pressure 
contribution. 

A problem similar to 1101) . l|llf) was studied in non- 
linear optics long ago pJJ for the opposite sign of the 
"pressure" p x in (|1U|) which corresponds to evolution of 
an optical beam in a focusing Kerr medium. More com- 
plex initial conditions were considered in [ill Il2j| . and 
in the recent paper [jjj the same nonlinear geometrical 
optics method has been applied to investigation of wave 
breaking phenomena in BEC. Here we shall apply the 
method of ^(j to the problem of BEC expansion. 

We look for a solution of ifTUjl. l(TT|) in the form 



p(x,t) 



v 2 x 2 



2f(t) v 2pf(ty 



, v(x,t) = x<j>(t), (12) 



where according to (|ll|l the functions fit) and 4>(t) must 
satisfy the initial conditions 



/(o) = i, m = o. 



(13) 



Substitution of (|12|l into (|10f) gives the relationship be- 
tween f(t) and 4>(t), 



2/(t)' 



(14) 
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as well as the differential equation for f(t) 

ff" = 2v 2 - (15) 

This equation can be readily solved with the initial con- 
ditions /(0) = 1, /'(0) = 2/(0)0(0) = (see JEJ and 
JHJ)) to give 



4i/t =2V/(/-l) 



In 



2/ - 1 + 2^/(7-1) 



(16) 



This formula determines implicitly / as a function of t 
and, hence, the distribution of density p(x,t) according 
to Eq. (|12|) . Then the function <j>(t) is determined by 
l|14jl. so that distribution of velocities v(x,t) = x<f>{t) is 
given by 



v(x, t) = 



/(*)' 



T- 



/(*)' 



(17) 



Thus, Eqs. (jP^)) . l(Ti|) . and (jTT))l give the complete analytic 
solution of the posed problem in the TF approximation. 

One of the effects which accompanies the expansion 
of the condensate is a non-monotonic behavior of the 
current density 



J{x,t) = p(x,t)v{x,t). 



(18) 



Equating its derivatives with respect to x and t to zero 
yields with the use of Eqs. (|TUj) . lfT2l and (|T?|) its maxi- 
mum value as well as corresponding values of x and t: 



J m («£m ) tm) 



(2p) 



3/2 



3/Lt 

Xrn \ I ~ T , t 

2v z 



27 ' 

\/3 + ln(2 + V3) 



(19) 



which with the use of Eq. JSJ can be expressed in terms 
of experimentally measurable parameters 



Jfn (*Em j tjn ) — V 

x m = 2.7 A v- 2 ! 3 



(20) 



Afa s ao 
~ S~ ' 

\ 1/3 
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V3 + ln(2 + V3)^ 76 ^ 
4^ 



The maximum of the density flow is proportional to 
the number of particles Af, although the time coordi- 
nate of the current maximum does not depend on the 
density but only on the condensate aspect ratio v. Let 
us consider the typical experiments on the condensate 
with M = 10 5 atoms of 87 Rb with scattering length 
a s ps 5nm 0] ■ Taking the length scale of the conden- 
sate of order of do « lpm (what corresponds to the fre- 
quency ujq ~ 5 ■ 10 3 Hz) and transversal radius ~ 10/im 
one can obtain for v = 0.2 the maximum of the current 
density J m ss 100 atoms-/im _2 -ms _1 in the coordinate 



x m w 11.5/xm (initial TF radius was ss 13.3/im) and at 
the time t rn w 1.1ms. 

In asymptotic limit of large t, when f(t) ^> /(0), we 
have 



f(t) = 2vt, (f>(t)^ 



1 

2?' 



Hence, the solution 1)1 2|l takes the form 

,•2 



p(x,t) = 



Avt 



1 
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v(x, t) = 



2f 



(21) 



(22) 



These formulas describe a hydrodynamic flow "by iner- 
tia" when the density becomes so small that the "pres- 
sure" does not accelerate the gas anymore (factor 1/2 in 
the second formula l|22l) corresponds to definition of the 
"time" variable in Eqs. (|TU|l ). 

From 12211 we can easily find asymptotic distribution 
of particles on their velocities: 



w{v) d v = l(i- v -y,. 



(23) 



which gives the number of atoms with velocities in the 
interval (v,v + dv). Also, in the framework of the TF ap- 
proximation Eq. (|12|) and with the use of (|15fl one obtains 
the velocity of the condensate front 



"/(*) 



p ld£ 

2 v dt 



2p\l- 



1/2 



(24) 



which at t — > oo asymptotically goes to vj(oo) = \/2p. 
The corresponding plots are shown in Fig. ^ 




FIG. 1: Time dependence of the front velocity Vf(t) for /i = 2 
and for different v (solid line v — 0.1, dotted line v — 0.2, 
dashed line v = 0.5). In the insert the corresponding initial 
distributions of the condensate density po(x) are depicted. 

To verify above findings, we have carried out numeri- 
cal simulations of the condensate dynamics governed by 
Eq. after switching off the trap potential, i.e. with 
v = 0. Initial profile of the condensate was taken as a 
numerical solution of Eq. (J5J . 
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FIG. 2: Spatial-temporal behavior of the density flow J(x,t) 
for n — 2 and v = 0.2. 

In agreement with analytical predictions we found 
that during expansion the current density J(x, t) = 
p(x,t)v(x,t) has maximum (see Fig. |2J). The space and 
time coordinates of this maximum which were calculated 
with the use of the TF approximation and by direct nu- 
merical calculations practically coincide (see Fig.|HJ). The 
discrepancy between the analytical and numerical calcu- 
lations increases with increase of the parameter v (see 
Fig. 3) what can be explained by lost of the accuracy of 
the TF approximation in accordance with the criterion 
M > v. 




X t 



FIG. 3: Dependence of the maximum of the density flow with 
coordinate (a) and time (b). Here solid thick lines correspond 
to the analytical solution and thin lines show the results of 
numerical calculations. 

Let us estimate time when we can consider an expan- 
tion of the condensate cloud only in the axial direction. 
In the particular case of the pancake geometry this can 



be done by conparison of the respect kinetic parts of 
the initial 3D GP equation. This gives us a criterium 
T /Ti « a 2 ja\ « 0.01 < 1 where T and T± are char- 
acteristics times of the processes in the axial and radial 
direction of the condensate correspondingly. Consider- 
ing time scale ( « Ti one can neglect the kinetics of 
the condensate in the radial direction and thus consider 
it to be unchanged in this direction. In order to esti- 
mate To and Tj_ we notice that velocity has an order of 
magnitude v ~ y/JI (see (24)) and thus T ~ ^fe, a± 
being the transverse radius. In the case of numerical cal- 
culations reported in Figs. 1,2 one estimates To ~ 0.4ms 
(1 dimensional unit) and T± ~ 40ms (100 dimmensional 
units) and thus the asymptotic limit described by (21) is 
indeed achieved. 

To conclude, it is interesting to make a comparison of 
our simplified ID model with some experiments where 
the expansion of the condensate cloud was observed. To 
follow the real experimental setup we take fi =const what 
corresponds to constant density po at the condensate cen- 
ter (see Eq. Jjj)) and look for the front velocities Vf(t) 
with different v what corresponds to different sizes (num- 
ber of particles AT) of the condensate in the transversal 
directions as was shown in Fig^ For larger values v 
(more narrow profile of the condensate) the front veloc- 
ity Vf is higher than for condensate with smaller v. This 
result of the ID model is in qualitative agreement with 
the experimental observations where faster expansion of 
the condensate cloud was observed in the direction with 
smaller transverse size 0, ■ 
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